In the interest of reproducibility, and to showcase our new package flotilla, I’ve reproduced many figures from the landmark single-cell paper, Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells by Shalek and Sujita, et al. Nature (2013).
Before we begin, let’s import everything we need.
In [1]:
# Import the flotilla package for biological data analysis
import flotilla
# Import "numerical python" library for number crunching
import numpy as np
# Import "panel data analysis" library for tabular data
import pandas as pd
# Import statistical data visualization package
# Note: As of November 6th, 2014, you will need the "master" version of
# seaborn on github (v0.5.dev), installed via
# "pip install git+ssh://git@github.com/mwaskom/seaborn.git
import seaborn as sns
# Turn on inline plots with IPython
%matplotlib inline
Shalek and Sujita, et al (2013)¶
In the 2013 paper, Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells (Shalek and Sujita, et al. Nature (2013)), Regev and colleagues performed single-cell sequencing 18 bone marrow-derived dendritic cells (BMDCs), in addition to 3 pooled samples.
Expression data¶
First, we will read in the expression data. These data were obtained using,
In [4]:
! wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE41nnn/GSE41265/suppl/GSE41265_allGenesTPM.txt.gz
We will also compare to the supplementary table 2 data, obtained using
In [5]:
! wget http://www.nature.com/nature/journal/v498/n7453/extref/nature12172-s1.zip
! unzip nature12172-s1.zip
In [6]:
expression = pd.read_table("GSE41265_allGenesTPM.txt.gz", compression="gzip", index_col=0)
expression.head()
Out[6]:
These data are in the “transcripts per million,” aka TPM unit. See this blog post if that sounds weird to you.
These data are formatted with samples on the columns, and genes on the rows. But we want the opposite, with samples on the rows and genes on the columns. This follows scikit-learn‘s standard of data matrices with size (n_samples, n_features) as each gene is a feature. So we will simply transpose this.
In [7]:
expression = expression.T
expression.head()
Out[7]:
The authors filtered the expression data based on having at least 3 single cells express genes with at TPM (transcripts per million, ) > 1. We can express this in using the pandas DataFrames easily.
First, from reading the paper and looking at the data, I know there are 18 single cells, and there are 18 samples that start with the letter “S.” So I will extract the single samples from the index (row names) using a lambda, a tiny function which in this case, tells me whether or not that sample id begins with the letter “S”.
In [8]:
singles_ids = expression.index[expression.index.map(lambda x: x.startswith('S'))]
print('number of single cells:', len(singles_ids))
singles = expression.ix[singles_ids]
expression_filtered = expression.ix[:, singles[singles > 1].count() >= 3]
expression_filtered = np.log(expression_filtered + 1)
expression_filtered.shape
Out[8]:
Hmm, that’s strange. The paper states that they had 6313 genes after filtering, but I get 6312. Even using “singles >= 1” doesn’t help.
(I also tried this with the expression table provided in the supplementary data as “SupplementaryTable2.xlsx,” and got the same results.)
Now that we’ve taken care of importing and filtering the expression data, let’s do the feature data of the expression data.
Expression feature data¶
This is similar to the fData from BioconductoR, where there’s some additional data on your features that you want to look at. They uploaded information about the features in their OTHER expression matrix, uploaded as a supplementary file, Supplementary_Table2.xlsx.
Notice that this is a csv and not an xlsx. This is because Excel mangled the gene IDS that started with 201* and assumed they were dates :(
The workaround I did was to add another column to the sheet with the formula ="'" & A1, press Command-Shift-End to select the end of the rows, and then do Ctrl-D to “fill down” to the bottom (thanks to this stackoverflow post for teaching me how to Excel). Then, I saved the file as a csv for maximum portability and compatibility.
In [10]:
expression2 = pd.read_csv('nature12172-s1/Supplementary_Table2.csv',
# Need to specify the index column as both the first and the last columns,
# Because the last column is the "Gene Category"
index_col=[0, -1], parse_dates=False, infer_datetime_format=False)
# This was also in features x samples format, so we need to transpose
expression2 = expression2.T
expression2.head()
Out[10]:
Now we need to strip the single-quote I added to all the gene names:
In [11]:
new_index, indexer = expression2.columns.reindex(map(lambda x: (x[0].lstrip("'"), x[1]), expression2.columns.values))
expression2.columns = new_index
expression2.head()
Out[11]:
We want to create a pandas.DataFrame from the “Gene Category” row for our expression_feature_data, which we will do via:
In [12]:
gene_ids, gene_category = zip(*expression2.columns.values)
gene_categories = pd.Series(gene_category, index=gene_ids, name='gene_category')
gene_categories
Out[12]:
In [13]:
expression_feature_data = pd.DataFrame(gene_categories)
expression_feature_data.head()
Out[13]:
Splicing Data¶
We obtain the splicing data from this study from the supplementary information, specifically the Supplementary_Table4.xls
In [14]:
splicing = pd.read_excel('nature12172-s1/Supplementary_Table4.xls', 'splicingTable.txt', index_col=(0,1))
splicing.head()
Out[14]:
In [15]:
splicing = splicing.T
splicing
Out[15]:
The three pooled samples aren’t named consistently with the expression data, so we have to fix that.
In [16]:
splicing.index[splicing.index.map(lambda x: 'P' in x)]
Out[16]:
Since the pooled sample IDs are inconsistent with the expression data, we have to change them. We can get the “P” and the number after that using regular expressions, called re in the Python standard library, e.g.:
In [17]:
import re
re.search(r'P\d', '10,000 cell Rep1 (P1)').group()
Out[17]:
In [18]:
def long_pooled_name_to_short(x):
if 'P' not in x:
return x
else:
return re.search(r'P\d', x).group()
splicing.index.map(long_pooled_name_to_short)
Out[18]:
And now we assign this new index as our index to the splicing dataframe
In [19]:
splicing.index = splicing.index.map(long_pooled_name_to_short)
splicing.head()
Out[19]:
Metadata
Now let’s get into creating a metadata dataframe. We’ll use the index from the expression_filtered data to create the minimum required column, 'phenotype', which has the name of the phenotype of that cell. And we’ll also add the column 'pooled' to indicate whether this sample is pooled or not.
In [20]:
metadata = pd.DataFrame(index=expression_filtered.index)
metadata['phenotype'] = 'BDMC'
metadata['pooled'] = metadata.index.map(lambda x: x.startswith('P'))
metadata
Out[20]:
Mapping stats data
In [21]:
mapping_stats = pd.read_excel('nature12172-s1/Supplementary_Table1.xls', sheetname='SuppTable1 2.txt')
mapping_stats
Out[21]:
Create a flotilla Study!
In [22]:
study = flotilla.Study(# The metadata describing phenotype and pooled samples
metadata,
# A version for this data
version='0.1.0',
# Dataframe of the filtered expression data
expression_data=expression_filtered,
# Dataframe of the feature data of the genes
expression_feature_data=expression_feature_data,
# Dataframe of the splicing data
splicing_data=splicing,
# Dataframe of the mapping stats data
mapping_stats_data=mapping_stats,
# Which column in "mapping_stats" has the number of reads
mapping_stats_number_mapped_col='PF_READS')
As a side note, you can save this study to disk now, so you can “embark” later:
In [23]:
study.save('shalek2013')
Note that this is saved to my home directory, in ~/flotilla_projects/<study_name>/. This will be saved in your home directory, too.
The datapackage.json file is what holds all the information relative to the study, and loosely follows the datapackage spec created by the Open Knowledge Foundation.
In [24]:
cat /Users/olga/flotilla_projects/shalek2013/datapackage.json
One thing to note is that when you save, the version number is bumped up. study.version (the one we just made) is 0.1.0, but the one we saved is 0.1.1, since we could have made some changes to the data.
Let’s look at what else is in this folder:
In [25]:
ls /Users/olga/flotilla_projects/shalek2013
So this is where all the other files are. Good to know!
We can “embark” on this newly-saved study now very painlessly, without having to open and process all those files again:
In [26]:
study2 = flotilla.embark('shalek2013')
Now we can start creating figures!
Figure 1
Here, we will attempt to re-create the sub-panels in Figure 1, where the original is:
Figure 1a
In [27]:
study.plot_two_samples('P1', 'P2')
Without flotilla, you would do
In [28]:
import seaborn as sns
sns.set_style('ticks')
x = expression_filtered.ix['P1']
y = expression_filtered.ix['P2']
jointgrid = sns.jointplot(x, y, joint_kws=dict(alpha=0.5))
xmin, xmax, ymin, ymax = jointgrid.ax_joint.axis()
jointgrid.ax_joint.set_xlim(0, xmax)
jointgrid.ax_joint.set_ylim(0, ymax)
Out[28]:
Figure 1b
Paper: $r=0.54$. Not sure at all what’s going on here.
In [29]:
study.plot_two_samples('S1', 'S2')
Without flotilla
In [30]:
import seaborn as sns
sns.set_style('ticks')
x = expression_filtered.ix['S1']
y = expression_filtered.ix['S2']
jointgrid = sns.jointplot(x, y, joint_kws=dict(alpha=0.5))
# Adjust xmin, ymin to 0
xmin, xmax, ymin, ymax = jointgrid.ax_joint.axis()
jointgrid.ax_joint.set_xlim(0, xmax)
jointgrid.ax_joint.set_ylim(0, ymax)
Out[30]:
By the way, you can do other kinds of plots with flotilla, like a kernel density estimate (“kde“) plot:
In [31]:
study.plot_two_samples('S1', 'S2', kind='kde')
Or a binned hexagon plot (“hexbin"):
In [32]:
study.plot_two_samples('S1', 'S2', kind='hexbin')
Any inputs that are valid to seaborn‘s jointplot are valid.
Figure 1c
In [33]:
x = study.expression.data.ix['P1']
y = study.expression.singles.mean()
y.name = "Average singles"
jointgrid = sns.jointplot(x, y, joint_kws=dict(alpha=0.5))
# Adjust xmin, ymin to 0
xmin, xmax, ymin, ymax = jointgrid.ax_joint.axis()
jointgrid.ax_joint.set_xlim(0, xmax)
jointgrid.ax_joint.set_ylim(0, ymax)
Out[33]:
Figure 2
Next, we will attempt to recreate the figures from Figure 2:
Figure 2a
For this figure, we will need the “LPS Response” and “Housekeeping” gene annotations, which weren’t very trivial to obtain, so I’ve moved them to the Appendix.
In [34]:
# Get colors for plotting the gene categories
dark2 = sns.color_palette('Dark2')
singles = study.expression.singles
# Get only gene categories for genes in the singles data
singles, gene_categories = singles.align(study.expression.feature_data.gene_category, join='left', axis=1)
mean = singles.mean()
std = singles.std()
jointgrid = sns.jointplot(mean, std, color='#262626', joint_kws=dict(alpha=0.5))
for i, (category, s) in enumerate(gene_categories.groupby(gene_categories)):
jointgrid.ax_joint.plot(mean[s.index], std[s.index], 'o', color=dark2[i], markersize=5)
jointgrid.ax_joint.set_xlabel('Standard deviation in single cells $\mu$')
jointgrid.ax_joint.set_ylabel('Average expression in single cells $\sigma$')
xmin, xmax, ymin, ymax = jointgrid.ax_joint.axis()
vmax = max(xmax, ymax)
vmin = min(xmin, ymin)
jointgrid.ax_joint.plot([0, vmax], [0, vmax], color='steelblue')
jointgrid.ax_joint.plot([0, vmax], [0, .25*vmax], color='grey')
jointgrid.ax_joint.set_xlim(0, xmax)
jointgrid.ax_joint.set_ylim(0, ymax)
jointgrid.ax_joint.fill_betweenx((ymin, ymax), 0, np.log(250), alpha=0.5, zorder=-1)
Out[34]:
I couldn’t find the data for the hESCs for the right-side panel of Fig. 2a, so I couldn’t remake the figure.
Figure 2b
In the paper, they use “522 most highly expressed genes (single-cell average TPM > 250)”, but I wasn’t able to replicate their numbers. If I use the pre-filtered expression data that I fed into flotilla, then I get 297 genes:
In [35]:
mean = study.expression.singles.mean()
highly_expressed_genes = mean.index[mean > np.log(250 + 1)]
len(highly_expressed_genes)
Out[35]:
Which is much less. If I use the original, unfiltered data, then I get the “522” number, but this seems strange because they did the filtering step of “discarded genes not appreciably expressed (transcripts per million (TPM) > 1) in at least three individual cells, retaining 6,313 genes for further analysis”, and yet they went back to the original data to get this new subset.
In [36]:
expression.ix[:, expression.ix[singles_ids].mean() > 250].shape
Out[36]:
In [37]:
expression_highly_expressed = np.log(expression.ix[singles_ids, expression.ix[singles_ids].mean() > 250] + 1)
mean = expression_highly_expressed.mean()
std = expression_highly_expressed.std()
mean_bins = pd.cut(mean, bins=np.arange(0, 11, 1))
# Coefficient of variation
cv = std/mean
cv.sort()
genes = mean.index
# for name, df in shalek2013.expression.singles.groupby(dict(zip(genes, mean_bins)), axis=1):
def calculate_cells_per_tpm_per_cv(df, cv):
df = df[df > 1]
df_aligned, cv_aligned = df.align(cv, join='inner', axis=1)
cv_aligned.sort()
n_cells = pd.Series(0, index=cv.index)
n_cells[cv_aligned.index] = df_aligned.ix[:, cv_aligned.index].count()
return n_cells
grouped = expression_highly_expressed.groupby(dict(zip(genes, mean_bins)), axis=1)
cells_per_tpm_per_cv = grouped.apply(calculate_cells_per_tpm_per_cv, cv=cv)
Here’s how you would make the original figure from the paper:
In [38]:
fig, ax = plt.subplots(figsize=(10, 10))
sns.heatmap(cells_per_tpm_per_cv, linewidth=0, ax=ax, yticklabels=False)
ax.set_yticks([])
ax.set_xlabel('ln(TPM, binned)')
Out[38]:
Doesn’t quite look the same. Maybe the y-axis labels were opposite, and higher up on the y-axis was less variant? Because I see a blob of color for (1,2] TPM (by the way, the figure in the paper is not TPM+1 as previous figures were)
This is how you would make a modified version of the figure, which also plots the coefficient of variation on a side-plot, which I like because it shows the CV changes directly on the heatmap. Also, technically this is $\ln$(TPM+1).
In [39]:
from matplotlib import gridspec
fig = plt.figure(figsize=(12, 10))
gs = gridspec.GridSpec(1, 2, wspace=0.01, hspace=0.01, width_ratios=[.2, 1])
cv_ax = fig.add_subplot(gs[0, 0])
heatmap_ax = fig.add_subplot(gs[0, 1])
sns.heatmap(cells_per_tpm_per_cv, linewidth=0, ax=heatmap_ax)
heatmap_ax.set_yticks([])
heatmap_ax.set_xlabel('$\ln$(TPM+1), binned')
y = np.arange(cv.shape[0])
cv_ax.set_xscale('log')
cv_ax.plot(cv, y, color='#262626')
cv_ax.fill_betweenx(cv, np.zeros(cv.shape), y, color='#262626', alpha=0.5)
cv_ax.set_ylim(0, y.max())
cv_ax.set_xlabel('CV = $\mu/\sigma$')
cv_ax.set_yticks([])
sns.despine(ax=cv_ax, left=True, right=False)
Figure 3
We will attempt to re-create the sub-panel figures from Figure 3:
Since we can’t re-do the microscopy (Figure 3a) or the RNA-FISH counts (Figure 3c), we will make Figures 3b. These histograms are simple to do outside of flotilla, so we do not have them within flotilla.
Figure 3b, top panel
In [40]:
fig, ax = plt.subplots()
sns.distplot(study.splicing.singles.values.flat, bins=np.arange(0, 1.05, 0.05), ax=ax)
ax.set_xlim(0, 1)
sns.despine()
Figure 3b, bottom panel
In [41]:
fig, ax = plt.subplots()
sns.distplot(study.splicing.pooled.values.flat, bins=np.arange(0, 1.05, 0.05), ax=ax, color='grey')
ax.set_xlim(0, 1)
sns.despine()
Figure 4a
Here, we can use the “interactive_pca” function we have to explore different dimensionality reductions on the data.
In [42]:
study.interactive_pca()
Out[42]:
A “sequences shortened” version of this is available as a gif:
Equivalently, I could have written out the plotting command by hand, instead of using study.interactive_pca:
In [43]:
study.plot_pca(feature_subset='gene_category: LPS Response', sample_subset='~pooled', plot_violins=False)
Out[43]:
Without flotilla, plot_pca is quite a bit of code:
In [44]:
import sys
from collections import defaultdict
from itertools import cycle
import math
from sklearn import decomposition
from sklearn.preprocessing import StandardScaler
import pandas as pd
from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from flotilla.visualize.color import dark2
from flotilla.visualize.generic import violinplot
class DataFrameReducerBase(object):
"""
Just like scikit-learn's reducers, but with prettied up DataFrames.
"""
def __init__(self, df, n_components=None, **decomposer_kwargs):
# This magically initializes the reducer like DataFramePCA or DataFrameNMF
if df.shape[1] <= 3:
raise ValueError(
"Too few features (n={}) to reduce".format(df.shape[1]))
super(DataFrameReducerBase, self).__init__(n_components=n_components,
**decomposer_kwargs)
self.reduced_space = self.fit_transform(df)
def relabel_pcs(self, x):
return "pc_" + str(int(x) + 1)
def fit(self, X):
try:
assert type(X) == pd.DataFrame
except AssertionError:
sys.stdout.write("Try again as a pandas DataFrame")
raise ValueError('Input X was not a pandas DataFrame, '
'was of type {} instead'.format(str(type(X))))
self.X = X
super(DataFrameReducerBase, self).fit(X)
self.components_ = pd.DataFrame(self.components_,
columns=self.X.columns).rename_axis(
self.relabel_pcs, 0)
try:
self.explained_variance_ = pd.Series(
self.explained_variance_).rename_axis(self.relabel_pcs, 0)
self.explained_variance_ratio_ = pd.Series(
self.explained_variance_ratio_).rename_axis(self.relabel_pcs,
0)
except AttributeError:
pass
return self
def transform(self, X):
component_space = super(DataFrameReducerBase, self).transform(X)
if type(self.X) == pd.DataFrame:
component_space = pd.DataFrame(component_space,
index=X.index).rename_axis(
self.relabel_pcs, 1)
return component_space
def fit_transform(self, X):
try:
assert type(X) == pd.DataFrame
except:
sys.stdout.write("Try again as a pandas DataFrame")
raise ValueError('Input X was not a pandas DataFrame, '
'was of type {} instead'.format(str(type(X))))
self.fit(X)
return self.transform(X)
class DataFramePCA(DataFrameReducerBase, decomposition.PCA):
pass
class DataFrameNMF(DataFrameReducerBase, decomposition.NMF):
def fit(self, X):
"""
duplicated fit code for DataFrameNMF because sklearn's NMF cheats for
efficiency and calls fit_transform. MRO resolves the closest
(in this package)
_single_fit_transform first and so there's a recursion error:
def fit(self, X, y=None, **params):
self._single_fit_transform(X, **params)
return self
"""
try:
assert type(X) == pd.DataFrame
except:
sys.stdout.write("Try again as a pandas DataFrame")
raise ValueError('Input X was not a pandas DataFrame, '
'was of type {} instead'.format(str(type(X))))
self.X = X
# notice this is fit_transform, not fit
super(decomposition.NMF, self).fit_transform(X)
self.components_ = pd.DataFrame(self.components_,
columns=self.X.columns).rename_axis(
self.relabel_pcs, 0)
return self
class DataFrameICA(DataFrameReducerBase, decomposition.FastICA):
pass
class DecompositionViz(object):
"""
Plots the reduced space from a decomposed dataset. Does not perform any
reductions of its own
"""
def __init__(self, reduced_space, components_,
explained_variance_ratio_,
feature_renamer=None, groupby=None,
singles=None, pooled=None, outliers=None,
featurewise=False,
order=None, violinplot_kws=None,
data_type='expression', label_to_color=None,
label_to_marker=None,
scale_by_variance=True, x_pc='pc_1',
y_pc='pc_2', n_vectors=20, distance='L1',
n_top_pc_features=50, max_char_width=30):
"""Plot the results of a decomposition visualization
Parameters
----------
reduced_space : pandas.DataFrame
A (n_samples, n_dimensions) DataFrame of the post-dimensionality
reduction data
components_ : pandas.DataFrame
A (n_features, n_dimensions) DataFrame of how much each feature
contributes to the components (trailing underscore to be
consistent with scikit-learn)
explained_variance_ratio_ : pandas.Series
A (n_dimensions,) Series of how much variance each component
explains. (trailing underscore to be consistent with scikit-learn)
feature_renamer : function, optional
A function which takes the name of the feature and renames it,
e.g. from an ENSEMBL ID to a HUGO known gene symbol. If not
provided, the original name is used.
groupby : mapping function | dict, optional
A mapping of the samples to a label, e.g. sample IDs to
phenotype, for the violinplots. If None, all samples are treated
the same and are colored the same.
singles : pandas.DataFrame, optional
For violinplots only. If provided and 'plot_violins' is True,
will plot the raw (not reduced) measurement values as violin plots.
pooled : pandas.DataFrame, optional
For violinplots only. If provided, pooled samples are plotted as
black dots within their label.
outliers : pandas.DataFrame, optional
For violinplots only. If provided, outlier samples are plotted as
a grey shadow within their label.
featurewise : bool, optional
If True, then the "samples" are features, e.g. genes instead of
samples, and the "features" are the samples, e.g. the cells
instead of the gene ids. Essentially, the transpose of the
original matrix. If True, then violins aren't plotted. (default
False)
order : list-like
The order of the labels for the violinplots, e.g. if the data is
from a differentiation timecourse, then this would be the labels
of the phenotypes, in the differentiation order.
violinplot_kws : dict
Any additional parameters to violinplot
data_type : 'expression' | 'splicing', optional
For violinplots only. The kind of data that was originally used
for the reduction. (default 'expression')
label_to_color : dict, optional
A mapping of the label, e.g. the phenotype, to the desired
plotting color (default None, auto-assigned with the groupby)
label_to_marker : dict, optional
A mapping of the label, e.g. the phenotype, to the desired
plotting symbol (default None, auto-assigned with the groupby)
scale_by_variance : bool, optional
If True, scale the x- and y-axes by their explained_variance_ratio_
(default True)
{x,y}_pc : str, optional
Principal component to plot on the x- and y-axis. (default "pc_1"
and "pc_2")
n_vectors : int, optional
Number of vectors to plot of the principal components. (default 20)
distance : 'L1' | 'L2', optional
The distance metric to use to plot the vector lengths. L1 is
"Cityblock", i.e. the sum of the x and y coordinates, and L2 is
the traditional Euclidean distance. (default "L1")
n_top_pc_features : int, optional
THe number of top features from the principal components to plot.
(default 50)
max_char_width : int, optional
Maximum character width of a feature name. Useful for crazy long
feature IDs like MISO IDs
"""
self.reduced_space = reduced_space
self.components_ = components_
self.explained_variance_ratio_ = explained_variance_ratio_
self.singles = singles
self.pooled = pooled
self.outliers = outliers
self.groupby = groupby
self.order = order
self.violinplot_kws = violinplot_kws if violinplot_kws is not None \
else {}
self.data_type = data_type
self.label_to_color = label_to_color
self.label_to_marker = label_to_marker
self.n_vectors = n_vectors
self.x_pc = x_pc
self.y_pc = y_pc
self.pcs = (self.x_pc, self.y_pc)
self.distance = distance
self.n_top_pc_features = n_top_pc_features
self.featurewise = featurewise
self.feature_renamer = feature_renamer
self.max_char_width = max_char_width
if self.label_to_color is None:
colors = cycle(dark2)
def color_factory():
return colors.next()
self.label_to_color = defaultdict(color_factory)
if self.label_to_marker is None:
markers = cycle(['o', '^', 's', 'v', '*', 'D', 'h'])
def marker_factory():
return markers.next()
self.label_to_marker = defaultdict(marker_factory)
if self.groupby is None:
self.groupby = dict.fromkeys(self.reduced_space.index, 'all')
self.grouped = self.reduced_space.groupby(self.groupby, axis=0)
if order is not None:
self.color_ordered = [self.label_to_color[x] for x in self.order]
else:
self.color_ordered = [self.label_to_color[x] for x in
self.grouped.groups]
self.loadings = self.components_.ix[[self.x_pc, self.y_pc]]
# Get the explained variance
if explained_variance_ratio_ is not None:
self.vars = explained_variance_ratio_[[self.x_pc, self.y_pc]]
else:
self.vars = pd.Series([1., 1.], index=[self.x_pc, self.y_pc])
if scale_by_variance:
self.loadings = self.loadings.multiply(self.vars, axis=0)
# sort features by magnitude/contribution to transformation
reduced_space = self.reduced_space[[self.x_pc, self.y_pc]]
farthest_sample = reduced_space.apply(np.linalg.norm, axis=0).max()
whole_space = self.loadings.apply(np.linalg.norm).max()
scale = .25 * farthest_sample / whole_space
self.loadings *= scale
ord = 2 if self.distance == 'L2' else 1
self.magnitudes = self.loadings.apply(np.linalg.norm, ord=ord)
self.magnitudes.sort(ascending=False)
self.top_features = set([])
self.pc_loadings_labels = {}
self.pc_loadings = {}
for pc in self.pcs:
x = self.components_.ix[pc].copy()
x.sort(ascending=True)
half_features = int(self.n_top_pc_features / 2)
if len(x) > self.n_top_pc_features:
a = x[:half_features]
b = x[-half_features:]
labels = np.r_[a.index, b.index]
self.pc_loadings[pc] = np.r_[a, b]
else:
labels = x.index
self.pc_loadings[pc] = x
self.pc_loadings_labels[pc] = labels
self.top_features.update(labels)
def __call__(self, ax=None, title='', plot_violins=True,
show_point_labels=False,
show_vectors=True,
show_vector_labels=True,
markersize=10, legend=True):
gs_x = 14
gs_y = 12
if ax is None:
self.reduced_fig, ax = plt.subplots(1, 1, figsize=(20, 10))
gs = GridSpec(gs_x, gs_y)
else:
gs = GridSpecFromSubplotSpec(gs_x, gs_y, ax.get_subplotspec())
self.reduced_fig = plt.gcf()
ax_components = plt.subplot(gs[:, :5])
ax_loading1 = plt.subplot(gs[:, 6:8])
ax_loading2 = plt.subplot(gs[:, 10:14])
self.plot_samples(show_point_labels=show_point_labels,
title=title, show_vectors=show_vectors,
show_vector_labels=show_vector_labels,
markersize=markersize, legend=legend,
ax=ax_components)
self.plot_loadings(pc=self.x_pc, ax=ax_loading1)
self.plot_loadings(pc=self.y_pc, ax=ax_loading2)
sns.despine()
self.reduced_fig.tight_layout()
if plot_violins and not self.featurewise and self.singles is not None:
self.plot_violins()
return self
def shorten(self, x):
if len(x) > self.max_char_width:
return '{}...'.format(x[:self.max_char_width])
else:
return x
def plot_samples(self, show_point_labels=True,
title='DataFramePCA', show_vectors=True,
show_vector_labels=True, markersize=10,
three_d=False, legend=True, ax=None):
"""
Given a pandas dataframe, performs DataFramePCA and plots the results in a
convenient single function.
Parameters
----------
groupby : groupby
How to group the samples by color/label
label_to_color : dict
Group labels to a matplotlib color E.g. if you've already chosen
specific colors to indicate a particular group. Otherwise will
auto-assign colors
label_to_marker : dict
Group labels to matplotlib marker
title : str
title of the plot
show_vectors : bool
Whether or not to draw the vectors indicating the supporting
principal components
show_vector_labels : bool
whether or not to draw the names of the vectors
show_point_labels : bool
Whether or not to label the scatter points
markersize : int
size of the scatter markers on the plot
text_group : list of str
Group names that you want labeled with text
three_d : bool
if you want hte plot in 3d (need to set up the axes beforehand)
Returns
-------
For each vector in data:
x, y, marker, distance
"""
if ax is None:
ax = plt.gca()
# Plot the samples
for name, df in self.grouped:
color = self.label_to_color[name]
marker = self.label_to_marker[name]
x = df[self.x_pc]
y = df[self.y_pc]
ax.plot(x, y, color=color, marker=marker, linestyle='None',
label=name, markersize=markersize, alpha=0.75,
markeredgewidth=.1)
try:
if not self.pooled.empty:
pooled_ids = x.index.intersection(self.pooled.index)
pooled_x, pooled_y = x[pooled_ids], y[pooled_ids]
ax.plot(pooled_x, pooled_y, 'o', color=color, marker=marker,
markeredgecolor='k', markeredgewidth=2,
label='{} pooled'.format(name),
markersize=markersize, alpha=0.75)
except AttributeError:
pass
try:
if not self.outliers.empty:
outlier_ids = x.index.intersection(self.outliers.index)
outlier_x, outlier_y = x[outlier_ids], y[outlier_ids]
ax.plot(outlier_x, outlier_y, 'o', color=color,
marker=marker,
markeredgecolor='lightgrey', markeredgewidth=5,
label='{} outlier'.format(name),
markersize=markersize, alpha=0.75)
except AttributeError:
pass
if show_point_labels:
for args in zip(x, y, df.index):
ax.text(*args)
# Plot vectors, if asked
if show_vectors:
for vector_label in self.magnitudes[:self.n_vectors].index:
x, y = self.loadings[vector_label]
ax.plot([0, x], [0, y], color='k', linewidth=1)
if show_vector_labels:
x_offset = math.copysign(5, x)
y_offset = math.copysign(5, y)
horizontalalignment = 'left' if x > 0 else 'right'
if self.feature_renamer is not None:
renamed = self.feature_renamer(vector_label)
else:
renamed = vector_label
ax.annotate(renamed, (x, y),
textcoords='offset points',
xytext=(x_offset, y_offset),
horizontalalignment=horizontalalignment)
# Label x and y axes
ax.set_xlabel(
'Principal Component {} (Explains {:.2f}% Of Variance)'.format(
str(self.x_pc), 100 * self.vars[self.x_pc]))
ax.set_ylabel(
'Principal Component {} (Explains {:.2f}% Of Variance)'.format(
str(self.y_pc), 100 * self.vars[self.y_pc]))
ax.set_title(title)
if legend:
ax.legend()
sns.despine()
def plot_loadings(self, pc='pc_1', n_features=50, ax=None):
loadings = self.pc_loadings[pc]
labels = self.pc_loadings_labels[pc]
if ax is None:
ax = plt.gca()
ax.plot(loadings, np.arange(loadings.shape[0]), 'o')
ax.set_yticks(np.arange(max(loadings.shape[0], n_features)))
ax.set_title("Component " + pc)
x_offset = max(loadings) * .05
ax.set_xlim(left=loadings.min() - x_offset,
right=loadings.max() + x_offset)
if self.feature_renamer is not None:
labels = map(self.feature_renamer, labels)
else:
labels = labels
labels = map(self.shorten, labels)
# ax.set_yticklabels(map(shorten, labels))
ax.set_yticklabels(labels)
for lab in ax.get_xticklabels():
lab.set_rotation(90)
sns.despine(ax=ax)
def plot_explained_variance(self, title="PCA explained variance"):
"""If the reducer is a form of PCA, then plot the explained variance
ratio by the components.
"""
# Plot the explained variance ratio
assert self.explained_variance_ratio_ is not None
import matplotlib.pyplot as plt
import seaborn as sns
fig, ax = plt.subplots()
ax.plot(self.explained_variance_ratio_, 'o-')
xticks = np.arange(len(self.explained_variance_ratio_))
ax.set_xticks(xticks)
ax.set_xticklabels(xticks + 1)
ax.set_xlabel('Principal component')
ax.set_ylabel('Fraction explained variance')
ax.set_title(title)
sns.despine()
def plot_violins(self):
"""Make violinplots of each feature
Must be called after plot_samples because it depends on the existence
of the "self.magnitudes" attribute.
"""
ncols = 4
nrows = 1
vector_labels = list(set(self.magnitudes[:self.n_vectors].index.union(
pd.Index(self.top_features))))
while ncols * nrows < len(vector_labels):
nrows += 1
self.violins_fig, axes = plt.subplots(nrows=nrows, ncols=ncols,
figsize=(4 * ncols, 4 * nrows))
if self.feature_renamer is not None:
renamed_vectors = map(self.feature_renamer, vector_labels)
else:
renamed_vectors = vector_labels
labels = [(y, x) for (y, x) in sorted(zip(renamed_vectors,
vector_labels))]
for (renamed, feature_id), ax in zip(labels, axes.flat):
singles = self.singles[feature_id] if self.singles is not None \
else None
pooled = self.pooled[feature_id] if self.pooled is not None else \
None
outliers = self.outliers[feature_id] if self.outliers is not None \
else None
title = '{}\n{}'.format(feature_id, renamed)
violinplot(singles, pooled_data=pooled, outliers=outliers,
groupby=self.groupby, color_ordered=self.color_ordered,
order=self.order, title=title,
ax=ax, data_type=self.data_type,
**self.violinplot_kws)
# Clear any unused axes
for ax in axes.flat:
# Check if the plotting space is empty
if len(ax.collections) == 0 or len(ax.lines) == 0:
ax.axis('off')
self.violins_fig.tight_layout()
# Notice we're using the original data, nothing from "study"
lps_response_genes = expression_feature_data.index[expression_feature_data.gene_category == 'LPS Response']
subset = expression_filtered.ix[singles_ids, lps_response_genes].dropna(how='all', axis=1)
subset_standardized = pd.DataFrame(StandardScaler().fit_transform(subset),
index=subset.index, columns=subset.columns)
pca = DataFramePCA(subset_standardized)
visualizer = DecompositionViz(pca.reduced_space, pca.components_, pca.explained_variance_ratio_)
visualizer()
Out[44]:
Figure 4b¶
In [45]:
lps_response_genes = study.expression.feature_subsets['gene_category: LPS Response']
lps_response = study.expression.singles.ix[:, lps_response_genes].dropna(how='all', axis=1)
lps_response.head()
Out[45]:
In [46]:
lps_response_corr = lps_response.corr()
“Elbow method” for determining number of clusters¶
The authors state that they used the “Elbow method” to determine the number of cluster centers. Essentially, you try a bunch of different $k$, and see where there is a flattening out of the metric, like an elbow. There’s a few different variations on which metric to use, such as using the average distance to the cluster center, or the explained variance. Let’s try the distance to cluster center first, because scikit-learn makes it easy.
In [48]:
from sklearn.cluster import KMeans
##### cluster data into K=1..10 clusters #####
ks = np.arange(1, 11).astype(int)
X = lps_response_corr.values
kmeans = [KMeans(n_clusters=k).fit(X) for k in ks]
# Scikit-learn makes this easy by computing the distance to the nearest center
dist_to_center = [km.inertia_ for km in kmeans]
fig, ax = plt.subplots()
ax.plot(ks, dist_to_center, 'o-')
ax.set_ylabel('Sum of distance to nearest cluster center')
sns.despine()
Not quite sure where the elbow is here. looks like there’s a big drop off after $k=1$, but that could just be an illusion. Since they didn’t specify which version of the elbow method they used, I’m not going to investigate this further, and just see if we can see what they see with the $k=5$ clusters that they found was optimal.
In [49]:
kmeans = KMeans(n_clusters=5)
lps_response_corr_clusters = kmeans.fit_predict(lps_response_corr.values)
lps_response_corr_clusters
Out[49]:
Now let’s create a dataframe with these genes in their cluster orders.
In [50]:
gene_to_cluster = dict(zip(lps_response_corr.columns, lps_response_corr_clusters))
dfs = []
for name, df in lps_response_corr.groupby(gene_to_cluster):
dfs.append(df)
lps_response_corr_ordered_by_clusters = pd.concat(dfs)
# Make symmetric, since we created this dataframe by smashing rows on top of each other, we need to reorder the columns
lps_response_corr_ordered_by_clusters = lps_response_corr_ordered_by_clusters.ix[:, lps_response_corr_ordered_by_clusters.index]
lps_response_corr_ordered_by_clusters.head()
Out[50]:
The next step is to get the principal-component reduced data, using only the LPS response genes. We can do this in flotilla using study.expression.reduce.
In [51]:
reduced = study.expression.reduce(singles_ids, feature_ids=lps_response_genes)
We can get the principal components using reduced.components_ (similar interface as scikit-learn).
In [52]:
reduced.components_.head()
Out[52]:
In [53]:
pc_components = reduced.components_.ix[:2, lps_response_corr_ordered_by_clusters.index].T
pc_components.head()
Out[53]:
In [54]:
import matplotlib as mpl
fig = plt.figure(figsize=(12, 10))
gs = gridspec.GridSpec(2, 2, wspace=0.1, hspace=0.1, width_ratios=[1, .2], height_ratios=[1, .1])
corr_ax = fig.add_subplot(gs[0, 0])
corr_cbar_ax = fig.add_subplot(gs[1, 0])
pc_ax = fig.add_subplot(gs[0, 1:])
pc_cbar_ax = fig.add_subplot(gs[1:, 1:])
sns.heatmap(lps_response_corr_ordered_by_clusters, linewidth=0, ax=corr_ax, cbar_ax=corr_cbar_ax,
cbar_kws=dict(orientation='horizontal'))
sns.heatmap(pc_components, cmap=mpl.cm.PRGn, linewidth=0, ax=pc_ax, cbar_ax=pc_cbar_ax,
cbar_kws=dict(orientation='horizontal'))
corr_ax.set_xlabel('')
corr_ax.set_ylabel('')
corr_ax.set_xticks([])
corr_ax.set_yticks([])
pc_ax.set_yticks([])
pc_ax.set_ylabel('')
Out[54]:
This looks pretty similar, maybe just rearranged cluster order. Let’s check what their data looks like when you plot this.
Their PC scores and clusters for the genes¶
In [55]:
gene_pc_clusters = pd.read_excel('nature12172-s1/Supplementary_Table5.xls', index_col=0)
gene_pc_clusters.head()
Out[55]:
In [56]:
data = lps_response_corr.ix[gene_pc_clusters.index, gene_pc_clusters.index].dropna(how='all', axis=0).dropna(how='all', axis=1)
fig = plt.figure(figsize=(12, 10))
gs = gridspec.GridSpec(2, 2, wspace=0.1, hspace=0.1, width_ratios=[1, .2], height_ratios=[1, .1])
corr_ax = fig.add_subplot(gs[0, 0])
corr_cbar_ax = fig.add_subplot(gs[1, 0])
pc_ax = fig.add_subplot(gs[0, 1:])
pc_cbar_ax = fig.add_subplot(gs[1:, 1:])
sns.heatmap(data, linewidth=0, square=True, vmin=-1, vmax=1, ax=corr_ax, cbar_ax=corr_cbar_ax, cbar_kws=dict(orientation='horizontal'))
sns.heatmap(gene_pc_clusters.ix[:, ['PC1 Score', 'PC2 Score']], linewidth=0, cmap=mpl.cm.PRGn,
ax=pc_ax, cbar_ax=pc_cbar_ax, cbar_kws=dict(orientation='horizontal'), xticklabels=False, yticklabels=False)
corr_ax.set_xlabel('')
corr_ax.set_ylabel('')
corr_ax.set_xticks([])
corr_ax.set_yticks([])
pc_ax.set_yticks([])
pc_ax.set_ylabel('')
Out[56]:
Sure enough, if I use their annotations, I get exactly that. Though there were two genes in their file that I didn’t have in the lps_response_corr data:
In [57]:
gene_pc_clusters.index.difference(lps_response_corr.index)
Oh joy, another datetime error, just like we had with expression2… Looking back at the original Excel file, there is one gene that Excel mangled to be a date:
Please, can we start using just plain ole .csvs for supplementary data! Excel does NOT preserve strings if they start with numbers, and instead thinks they are dates.
In [58]:
import collections
collections.Counter(gene_pc_clusters.index.map(type))
Out[58]:
Yep, it’s just that one that got mangled…. oh well.
In [59]:
gene_pc_clusters_genes = set(filter(lambda x: isinstance(x, unicode), gene_pc_clusters.index))
gene_pc_clusters_genes.difference(lps_response_corr.index)
Out[59]:
So, “RPS6KA2” is the only gene that was in their list of genes and not in mine.
Supplementary figures¶
Now we get to have even more fun by plotting the Supplementary figures! :D
Ironically, the supplementary figures are usually way easier to access (like not behind a paywall), and yet they’re usually the documents that really have the crucial information about the experiments.
Supplementary Figure 1¶
In [60]:
singles_mean = study.expression.singles.mean()
singles_mean.name = 'Single cell average'
# Need to convert "average_singles" to a DataFrame instead of a single-row Series
singles_mean = pd.DataFrame(singles_mean)
singles_mean.head()
Out[60]:
In [61]:
data_for_correlations = pd.concat([study.expression.singles, singles_mean.T, study.expression.pooled])
# Take the transpose of the data, because the plotting algorithm calculates correlations between columns,
# And we want the correlations between samples, not features
data_for_correlations = data_for_correlations.T
data_for_correlations.head()
# %time sns.corrplot(data_for_correlations)
Out[61]:
In [62]:
fig, ax = plt.subplots(figsize=(10, 10))
sns.corrplot(data_for_correlations, ax=ax)
sns.despine()
Notice that this is mostly red, while in the figure from the paper, it was both blue and red. This is because the colormap started at 0.2 (not negative), and was centered with white at about 0.6. I see that they’re trying to emphasize how much more correlated the pooled samples are to each other, but I think a simple sequential map would have been more effective.
Supplementary Figures 2 and 3¶
Supplementary Figure 4¶
Supplementary Figure 4 was from published data, however the citation in the Supplementary Information (#23) was a machine-learning book, and #23 in the main text citations was a review of probabilistic graphical models, neither of which have the mouse embryonic stem cells or mouse embryonic fibroblasts used in the figure.
Supplementary Figure 5¶
For this figure, we can only plot 5d, since it’s derived directly from a table in their dataset.
Warning: these data are going to require some serious cleaning. Yay data janitorial duties!
Supplementary Figure 5d¶
In [63]:
barcoded = pd.read_excel('nature12172-s1/Supplementary_Table7.xlsx')
barcoded.head()
Out[63]:
The first three columns are TPM calculated from the three samples that have molecular barcodes, and the last three columns are the integer counts of molecular barcodes from the three molecular barcode samples.
Let’s remove the “Unnamed: 3” column which is all NaNs. We’ll do that with the .dropna method, specifying axis=1 for columns and how="all" to make sure only columns that have ALL NaNs are removed.
In [64]:
barcoded = barcoded.dropna(how='all', axis=1)
barcoded.head()
Out[64]:
Next, let’s drop that pesky “GENE” row. Don’t worry, we’ll get the sample ID names back next.
In [65]:
barcoded = barcoded.drop('GENE', axis=0)
barcoded.head()
Out[65]:
We’ll create a pandas.MultiIndex from the tuples of (sample_id, measurement_type) pair.
In [66]:
columns = pd.MultiIndex.from_tuples([('MB_S1', 'TPM'),
('MB_S2', 'TPM'),
('MB_S3', 'TPM'),
('MB_S1', 'Unique Barcodes'),
('MB_S2', 'Unique Barcodes'),
('MB_S3', 'Unique Barcodes')])
barcoded.columns = columns
barcoded = barcoded.sort_index(axis=1)
barcoded.head()
Out[66]:
For the next move, we’re going to do some crazy pandas-fu. First we’re going to transpose, then reset_index of the transpose. Just so you know what this looks like, it’s this.
In [67]:
barcoded.T.reset_index().head()
Out[67]:
Next, we’re going to transform the data into a tidy format, with separate columns for sample ids, measurement types, the gene that was measured, and its measurement value.
In [68]:
barcoded_tidy = pd.melt(barcoded.T.reset_index(), id_vars=['level_0', 'level_1'])
barcoded_tidy.head()
Out[68]:
Now let’s rename these columns into something more useful, instead of “level_0”
In [69]:
barcoded_tidy = barcoded_tidy.rename(columns={'level_0': 'sample_id', 'level_1': 'measurement', 'variable': 'gene_name'})
barcoded_tidy.head()
Out[69]:
Next, we’re going to take some seemingly-duplicating steps, but trust me, it’ll make the data easier.
In [70]:
barcoded_tidy['TPM'] = barcoded_tidy.value[barcoded_tidy.measurement == 'TPM']
barcoded_tidy['Unique Barcodes'] = barcoded_tidy.value[barcoded_tidy.measurement == 'Unique Barcodes']
Fill the values of the “TPM“‘s forwards, since they appear first, and fill the values of the “Unique Barcodes” backwards, since they’re second
In [71]:
barcoded_tidy.TPM = barcoded_tidy.TPM.ffill()
barcoded_tidy['Unique Barcodes'] = barcoded_tidy['Unique Barcodes'].bfill()
barcoded_tidy.head()
Out[71]:
Drop the “measurement” column and drop duplicate rows.
In [72]:
barcoded_tidy = barcoded_tidy.drop('measurement', axis=1)
barcoded_tidy = barcoded_tidy.drop_duplicates()
barcoded_tidy.head()
Out[72]:
In [73]:
barcoded_tidy['log TPM'] = np.log(barcoded_tidy.TPM)
barcoded_tidy['log Unique Barcodes'] = np.log(barcoded_tidy['Unique Barcodes'])
Now we can use the convenient linear model plot (lmplot) in seaborn to plot these three samples together!
In [74]:
sns.lmplot('log TPM', 'log Unique Barcodes', barcoded_tidy, col='sample_id')
Out[74]:
Supplementary Figures 6-20¶
Conclusions¶
While there may be minor, undocumented, differences between the methods presented in the manuscript and the figures, the application of flotilla presents an opportunity to avoid these types of inconsistencies by strictly documenting every change to code and every transformation of the data. The biology the authors found is clearly real, as they did the knockout experiment of Ifnr-/- and saw that indeed the maturation process was affected, and Stat2 and Irf7 had much lower expression, as with the “maturing” cells in the data.
comments powered by Disqus